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We investigate the structure and equilibrium linear-response dynamics of suspensions of hard 
colloidal dumbbells using Brownian Dynamics computer simulations. The focus lies on the 
dense fluid and plastic crystal states of the colloids with investigated aspect (elongation-to- 
diameter) ratios varying from the hard sphere limit up to 0.39, which is roughly the stability 
limit of the plastic crystal phase. We find expected structural changes with larger elongation 
with respect to the hard sphere reference case and very localized orientational correlations, 
typically just involving next-neighbor couplings. These relatively weak correlations are also 
reflected in only minor effects on the translational and rotational diffusion coefficients for 
most of the investigated elongations. However, the linear response shear viscosity exhibits a 
dramatic increase at high packing fractions [(j) > 0.5) beyond a critical anisotropy factor of 
about L* ~ 0.15 which is surprising in view of the relatively weak changes found before on the 
level of colloidal self-dynamics. We suspect that even for the small investigated anisotropies, 
newly occurring, collective rotational-translational couplings must be made responsible for 
the slow time scales appearing in the plastic crystal. 

Keywords: dense fluids; plastic crystals; colloidal dumbbells; anisotropic colloids; linear 
response dynamics; Brownian dynamics simulation; 


1. Introduction 

Suspensions of hard spherical colloids have been extensively studied theoretically 
as well as experimentally as a model system for simple liquids in and out-of equilib¬ 
rium m- However, even weak anisotropies of the particle shape can have dramatic 
effects on the structure, dynamics, and phase behavior of colloidal suspensions 
while the details of their molecular origins are not well understood. Particularly 
interesting is the question at which elongation what dynamic and structural prop¬ 
erties start to deviate from the limiting hard sphere reference case. To address 
such a fundamental question, well-defined experimental and theoretical model sys¬ 
tems are in need. One of the simplest realizations of mildly anisotropic colloids are 
dumbbell-shaped particles (that is, a dimer of two fused colloidal spheres) which can 
be synthesized nowadays by different routes in well-controlled and monodisperse 
ways [THIS]. The phase diagram of hard-core dumbbells is theoretically known, 
where, at high packing fractions, the colloids can form the so-called plastic crystal 
phase |191425j . Here, the translational degrees of freedom are essentially frozen at 
high packings as in the hard sphere system but the particles are free to rotate. 
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Experimental realizations of dumbbell suspensions have been shown to fit well into 
the respective regions of the theoretical phase diagram [IHlEe]. Thanks to these 
recent advances, viable, purely repulsive monodisperse model systems are exper¬ 
imentally available which crystallize in equilibrium conditions without the help 
of external fields or long-range interactions. Colloidal dumbbells are thus among 
the most promising experimental model systems for a fundamental understanding 
of structural and dynamic effects of increasing the elongation with respect to the 
hard sphere reference case, in particular in the (plastic) crystal phase where new 
phenomena are expected m- 

In this work, we systematically study the structure and fluctuations of suspen¬ 
sions of steeply repulsive dumbbell particles in equilibrium by means of Brownian 
Dynamics (BD) computer simulations. We focus on the dense regime of the fluid (F) 
to plastic crystal (PC) transition region of the phase diagram which is well accessi¬ 
ble for experimental dumbbell suspensions and promises interesting insight into the 
influence of elongation on steric correlations and structural properties in both fluid 
and crystal. Apart from simple spatial correlation functions, we also systematically 
explore the dynamic properties in the linear response regime at high volume frac¬ 
tions measured by equilibrium fluctuations. All results for the elongated dumbbells 
are discussed with respect to the well-explored hard sphere reference system. 


2. Theoretical methods 


2.1. Dumbbell model and Brownian Dynamics 

In our model a dumbbell particle consists of two spherical beads of diameter D 
rigidly constrained at a center-to-center distance L, see Fig. With that we define 
a dimensionless elongation, or, aspect ratio, as L* = L/D. For L* = 0, we recover 
the reference case of one spherical colloid. 



Figure 1. Sketch of the dumbbell geometry with dimensionless aspect ratio L* = LfD. 

The site-site (bead-bead) pair interaction is modeled by a Yukawa (or ’screened 
Coulomb’) potential which interpolates between the long-ranged electrostatic 
Coulomb and the purely excluded volume interaction. The potential between two 
beads at distance r is given by 


V (r) = e— exp {—K (r — cj)} . (1) 

r 

The parameter k is a screening constant, which may be used to tune the softness 
of the pair interaction. In the limit of k —oo one obtains the purely hard-core 
interaction. We choose k = 20/cr as a compromise between a steep hard-core like 
repulsion at the same time being smooth enough to not cause numerical problems 
in the integration. The parameters a and e = ksT define the length and energy 
scale of the interaction in our model, respectively. Since we compare our results to 
the hard-dumbbell and hard-sphere reference system, we define the effective bead 
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diameter D via the Barker-Henderson relation m 


oo 



0 


( 2 ) 


The effective volume fraction (j) = p{7r/6)a{L*)D^ at dumbbell number density p is 
defined accordingly. The function a{L*) = 1 + 3L*/2 — L*^j2 is a geometric factor 
accounting for the bead overlap volume [20] . 

The total instantaneous force F ij (t) exerted by particle j on particle i is the sum 
of the forces F(j acting between the k-th site of the j-th particle and the 

l-th. site of the j-th particle, via 


2 

Fij(t) = F(j (t). (3) 

k,l=l 

The center of mass (COM) vector of the Tth particle is decomposed with respect 
to the axial symmetry as 


Ri(t) = Ri,_L(t) + Ri,||(t), (4) 

Rj_ll(t) = [rii(t) • Ri(t)] (5) 

where is the unit vector pointing in the direction of the connection between 

the centers of the beads belonging to the i-th particle. The total force Fj(t) is 
decomposed analogously: The parallel and perpendicular parts are given by 


Fi(t) = Fi^_L(t)+Fj_||(t), (6) 

Fi,||(t) = [a(t).R(t)]a(t). (7) 

Given this framework, we implement the Brownian dynamics (BD) algorithm 
for the above defined coordinates as follows: The particle positions are updated 
using an explicit forward Euler method with the finite time step At <C t. Start¬ 
ing from the configuration {Rr,nn at time t = nAt, we find new coordinates 
|Rn+i^ a,t t -|- At following the scheme 


R:^+i = R(;il+AtiF-||+5r,,|,nr, 


R”t^ — Ri(_L + 


( 8 ) 

(9) 


where and eip are unit vectors oriented perpendicular to the director fij. The 
random variates dvi^a, « £ {11)1)2} have zero mean and the variances = 

2T)|jAt and = 2D±At, respectively. 

The torque acting on symmetric dumbbells is given by Tj(t) = x 

(F(j 2 ),±(t) — F(j_i)^_i_(t)), where F(j fc) _i_(t) is the total instantaneous force on the 
fc-th bead of the z-th particle. The internal director of the z-th particle is updated 
according to 


fl 


n+1 _ 


= nr 




( 10 ) 
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where 5xi and 8 x 2 are zero mean Gaussian distributed random variates with vari¬ 
ance (^ 8 x'j'^ = 2DrAt. The single particle mobility is given by parallel Dy, per¬ 
pendicular D± and rotational Dj. diffusion coefficients. We neglect hydrodynamic 
interactions between the colloids which is justified for the high packing fractions 
where steric correlations dominate the equilibrium structure and linear response 
behavior [28l [29j . 

The infinite-dilution (single-particle) self-diffusion coefficients of a dumbbell par¬ 
ticle were calculated previously using the shell-bead model (SHM) [30H3l] . In this 
method, the particle surface is represented by a number of mini-beads which act as 
sources of hydrodynamic friction, where no-slip boundary conditions are assumed. 
The mobility tensors are calculated on the Oseen-level using an increasing number 
of mini-beads. Using these tensors, the values are extrapolated to zero mini-bead 
size, e.g., infinite number of friction sources on the particle surface. We are us¬ 
ing the parallel T>||/T>q, perpendicular D±/Dq and rotational Dr/D^ diffusion 
coefficients obtained by the above method as input for our BD simulations. The 
free single-sphere {L* = 0) translational Dq and rotational diffusivities may 
be obtained from the respective Stokes-Einstein relations. The translational COM 
diffusion coefficient of a free dumbbell is given by Dq = ^Z)|| -|- The Brownian 
time is accordingly given by r = <j^/Dq j35l [36] . 

In our simulations, the systems are subject to periodic boundary conditions in 
all three Cartesian dimensions. The presented results have been obtained from 
BD simulation runs with N = 864 particles in a constant volume V, so that the 
number density p = N/V. The total simulation time is lOOr with the time step of 
At = 10“^r. The face-centered cubic (fee) ordered COMs with directors pointing in 
the (1,1,1) directions have been set as initial configuration. To ensure equilibrium 
conditions results have only been obtained for t > 50t. All statistical correlations 
have been calculated for 50t < t < lOOr using every 100-th time step. 



Figure 2. Schematic representation of the phase diagram of hard dumbbells m- We investigate systems 
in the fluid (F) and the plastic crystal (PC) phase. The filled diamonds (♦) indicate the state points 
investigated in this work, according to a Barker-Henderson mapping § of the site-site Yukawa pair 
potential 0 onto hard spheres. 




July 10, 2015 


Molecular Physics dumbbells'final 


2.2. Static structure 

The expansion of the full radial distribution function |37j yields two partial corre¬ 
lation functions of particular interest: We calculate the radial distribution function 
of the centers of masses (COMs) 


\ i^j / 

and the orientational radial distribution function 

^ pN^ir) ^ (r - , (12) 

where Rjj denotes the COM separation vector of the i-th and j-th particles, and 
Oij is the angle between the respective directors. The orientational distribution 
function gp^{r) describes the spatial correlation of the particle directors with P 2 
denoting the second Legendre polynomial. The function gp^ir) equals 1 for parallel 
aligned particles, — ^ for perpendicular alignment, and vanishes for uncorrelated 
orientations. 


2.3. Dynamical properties 
2.3.1. Translational diffusion 

We calculate the translational diffusion in the linear response regime from Green- 
Kubo-type of relations. In the limit of linear response, we thus obtain the long-time 
COM self-diffusion coefficient from the time correlation functions in equilibrium 
from the integral of the velocity autocorrelation function (VACF) as 

t 

D{t) = I (V(0)-V(t'))dt', (13) 

0 

where V (t) is the instantaneous COM velocity. The long-time limit of this function 
is the long-time self-diffusion coefficient 

= lim D{t). (14) 


2.3.2. Orientational relaxation 

The linear response orientational relaxation is quantified in terms of the 
directional autocorrelation function (DACF) defined as 

Cn{t) = (fi(0) • n{t)) = (15) 

We assume that the particle orientations become diffusive for long times with the 
constant therefore asymptotically the following relation holds [36] : 

Cn{t) = e-^^^K 


(16) 
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2.3.3. Shear viscosity 

In the limit of zero strain, we may obtain the dynamic viscosity of the suspen¬ 
sion from equilibrium fluctuations. In particular, we calculate the relative dynamic 
viscosity difference from the off-diagonal stress autocorrelation function (SACF) 


where the symmetric stress tensor is dehned as 

^ N N 

~ Y f'ijaFijP, 

i=l j>i 


(17) 


(18) 


and a ^ (5 denote Cartesian components. From our BD simulations, the potential 


part of the stress tensor is available, therefore (18) does not contain a momentum 
part. 

The steady shear viscosity difference is obtained from the integral of the SACF 
as 




OO 

j Zal3it)dt. 

0 


(19) 


3. Results and discussion 
3.1. Spatial structure 

In Fig.|^the COM radial distribution functions (RDFs) of dumbbell suspensions at 
the volume fraction (j) = 0.60 are displayed for various elongations. The signature 
of the phase transition from PC to fluid (F) may be readily observed here. At 
state points with L* < 0.3 the system is in the PC phase, cf. Fig. where the 
COM positions are frozen in a crystal lattice and the directors do not show any 
long-distance correlations. Hence, the corresponding RDFs show strong correlations 
and long-ranged oscillations. They vanish at L* = 0.39, where the system is in the 
isotropic fluid phase. The inset of Fig. shows the center of mass (COM) radial 
distribution functions (RDFs) g{r) at the fixed colloidal packing fraction (p = 0.44. 
All shown curves are in the fluid (F) phase. With increasing elongation, the extrema 
get less pronounced and are shifted to greater distances due to the larger effective 
particle sizes. 

The orientational pair correlation functions (PCFs) gp^ at the volume fractions 
4> = 0.60 and p = 0.44 are depicted in Fig. and its inset, respectively, for various 
representative elongations. As expected, the orientational correlation in space is 
essentially flat for suspensions of nearly spherical particles {L* = 0.02) at all volume 
fractions. At roughly rja ~ 1, the gp.^ for non-zero elongations show maxima 
indicating a strong correlation of nearest neighbors at contact. This first correlation 
peak moves slightly away from r = cr as the elongation increases. At the hrst 
maxima, the functions are positive, representing a preferably parallel orientation 
of dumbbell particles at close contact. For a bit larger distances, r/cr ~ 1.1 — 1.3, 
negative dips indicating orientational anti-correlations are observed which tend to 
higher distances as the elongation is increased. Hence, in the PC phase at p = 0.60, 
there are non-vanishing next neighbor correlations which grow with elongation. 
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Figure 3. The radial distribution function g(r) for the dumbbell center-of-mass at a volume fraction 
4> = 0.60: black, blue, and purple are in the PC, and red in the F phases, respectively. Inset: The same for 
a packing fraction = 0.44, where all systems are in the F phase. 


With further increasing aspect ratio, the orientational correlation becomes non-zero 
over a higher distance as the system crosses to the dense fluid phase at 0 = 0.44. 
Here, we observe a new correlation peak at in the dense fluid state for L* = 
0.39 close to r/a ~ 1.8 and 2.8, indicating growing second and third neighbor 
correlations. Thus, in the fluid phase the correlations are longer ranged due to the 
higher disorder and collisions when compared to the PC phase. 



Figure 4. Orientational pair correlation function gp 2 {'f') at a volume fraction ^ = 0.60: the black, blue, 
and purple curves present data for the systems in the plastic crystal (PC) phase and red in the fluid (F) 
phase, respectively. Inset: The same for a packing fraction (p = 0.44 where all systems are in the fluid 
phase. 


3.2. Translational diffusion 

3.2.1. Time-dependent diffusion coefficient 

Figure a) and b) show the time-dependent COM diffusion coefficients D{t) 
at constant elongation L* = 0.19 and L* = 0.39 for different colloidal volume 
fractions, respectively. At early times {t < 10“^), before the particles feel the inter¬ 
acting neighbors in a non-dilute suspension, the diffusion coefficients are close to 




























July 10, 2015 


Molecular Physics dumbbells'final 



Figure 5. Time-dependent diffusion coefficient Z)(t)/Do from the VACF at (a) L* = 0.19 and (b) L* = 
0.39 for different packing fractions. The dashed line in (b) is at state point 4> = 0.44, L* = 0.02 for a direct 
comparison. 


the short-time limits of a single, free particle. Within times t < 10“^r, a continu¬ 
ous transition to long-time diffusion starts which is reached at about the Brownian 
time scale r. With increasing volume fraction </>, the cross-over from short-time 
to long-time behavior sets in earlier and is steeper. At the elongation L* = 0.19 
the systems with packing fractions </> = 0.54 and (p = 0.60 are in the PC phase. 
Due to the translational constraints in the crystal, the diffusion vanishes at about 
a tenth of a Brownian time. The comparison of the diffusion between L* = 0.19 
and L* = 0.39 at cp = 0.44 in the dense fluid in panel (a) demonstrates that the 
influence of the aspect ratio is very small in the range of focus of this study. Since 
all the shown curves for the elongation L* = 0.39 are in the F phase (panel (b)), 
the diffusion coefficients are not vanishing for all volume fractions. 

In Fig. 1^ the values of the long-time COM self-diffusion coefficients are pre¬ 
sented for various elongations and packing fractions. As expected, the diffusion is 
significantly decreased for larger packer fractions, as known for the spherical refer¬ 
ence case. Surprisingly, the dependence of the diffusion on elongation is weak and 
we find that this is the case for both parallel and perpendicular long-time diffusion 
as well. Apparently, the anisotropy is not large enough to substantially change the 
long-time (single colloid) friction in these systems. Furthermore, the data shows a 
discontinuous transition to vanishing diffusion at elongations where the transition 
to the PC phase takes place, see the curves for the large packing fractions cp > 0.54. 
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Figure 6. Long-time self-diffusion coefficients Df' obtained from the VACFs via relation 1^. 

3.3. Rotational relaxation and diffusion 

Fig. [7] shows the rotational diffusion coefficients obtained from exponential fits to 
the DACFs. At small elongations L* < 0.2 the rotational diffusion is almost con¬ 
stant with respect to volume fraction, albeit the system crosses from the fluid to 
the PC phase at high fractions. Hence, for these small elongations the particles 
indeed rotate almost freely in the PC phase. The obvious statistical outliers in the 
rotational diffusion data for small packing fractions (p < 0.3 are within the statis¬ 
tical error of ±5%, estimated from block averages over independent trajectories of 
length of lOr. However, suspensions of dumbbells with L* = 0.29 show a notable 
drop of the rotational diffusion coefficients for densities higher than p = 0.5 and 
packing effects clearly influence the rotation correlations in the PC phase. On in¬ 
creasing the elongation further, now a clear non-linear density-dependence of the 
rotational relaxation emerges which becomes obvious for L* = 0.39 for packings 
larger than p > 0.3. Here, the system is in the fluid phase only and the closer 
contacts in the disordered systems alter the rotational dynamics substantially. At 
the highest investigated packing fraction [p = 0.60) the rotational diffusion drops 
down by more than 30 %. 



Figure 7. Dependence of rotational diffusion coefficients Dp at various elongations from exponential fits 
to DACF on volume fraction 0. The statistical error of the data is estimated to be ±5 %, calculated from 
block averages over independent trajectories of length of lOr. 
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3.4. Shear viscosity 

Figure shows the stress autocorrelation function SACFs at various aspect ratios 
L*. For the smaller packing fractions cj) < 0.5, the correlation functions decay slower 
for higher densities while the elongation dependence is rather weak. However, at the 
volume fraction cp = 0.60 the system is in the PC phase for the shown elongations, 
except for L* = 0.39. Here an interesting behavior can be observed: at L* = 0.02 
the SACF decays more rapidly in the PC state point than in the fluid states. 
However, for the next two larger elongations L* = 0.19 and L* = 0.29, the SACF 
decays slower and develops a new slow time scale for times about 10~^r in the PC 
phase. This must be assigned to slow stress relaxations in the dense crystal phase 
likely stemming from more pronounced translational-rotational couplings m , since 
they are absent for the small elongations at the same density. This long-time tail 
occurring in the PC phase seems to diminish when crossing into the F phase by 
increasing the aspect ratio to L* = 0.39. 



Figure 8. SACFs at different elongations (increasing from top left to bottom right panels) and volume 
fractions. 


Figure a) shows the packing fraction dependence of the relative steady shear 
viscosity r/r-o = Vo/ Vs for different elongations L*, with r]s being the viscosity of the 
suspending medium. These data are obtained from integration of the SACFs in 
Fig. As an important comparison, also the hard sphere (HS) reference data fol¬ 
lowing the empirical Krieger-Dougherty relation [281 [38] for the purely fluid phase 
is shown. As we readily see, there is a difference of the viscosity for a fluid of hard 
spheres when compared to our minimal elongation data (at L* = 0.02). This can be 
understood as a result of the sensitivity of viscosities to the softness of the chosen 
pair potential for repulsive spheres [39]. In the present case, the L* = 0.02 system 
with minute anisotropy can be considered equal to the purely spherical system. 
At higher packing the HS system is in the crystal phase and the linear shear re¬ 
sponse drops to zero (and not accounted for in the fluid-state Krieger-Dougherty 
approach). The shear viscosity of dumbbells is similar to the HS system in the 
fluid phase below packings of < 0.45 but deviates strongly for larger packings 
for all elongations where the dumbbell shear response substantially increases with 
L*. We explain the increase by the long time tails in the SACFs in Figure]^ due 
to the relaxation of the rotational degrees of freedom missing in the HS case. The 
large increase of the shear response is somewhat unexpected since the translational 
degrees of freedom are still frozen, except for the largest elongation L* = 0.39 for 
which all shown data are in the fluid phase. We note that the error bars of this data 
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are hard to estimate due to the long-time tails in the SACFs. While we believe from 
the systematic behavior of the data that the trends are reliable, small fluctuations 
of the data as for cf) = 0.44 are very probable within the statistical uncertainty. 




Figure 9. Steady shear relative viscosity difference rjro (a) as a function of the volume fraction cp and 
(b) versus elongation L* at high densities. The hard sphere (HS) line ( ) refers to the phenomenological 
Krieger-Dougherty theory |281138| . 


Figure b) shows the dependence of the relative steady shear viscosity on the 
elongation at the high packing fractions (p > 0.39. This view reveals a transition 
from the hard spherical case of almost vanishing viscosity to a finite level at about 
a critical elongation of about L* = 0.15. This is most evident at the highest pack¬ 
ing fraction where the transition is relatively steep. (The fluctuations of the data 
between L* = 0.2 and 0.3 are most certainly within the statistical uncertainty.) 
Hence, it may be inferred from the equilibrium linear response data that a small 
deviation from the spherical shape already to an aspect ratio of L* > 0.15 induces 
a dramatic collective response that manifests itself in a large increase of the linear 
response shear viscosity. 


4. Summary and concluding remarks 

In summary, we have investigated the equilibrium structure and fluctuations of 
colloidal dumbbells by means of Brownian Dynamics computer simulations. We 
focused on high packing fractions and on the weak elongation regime {L* < 0.4) 
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where the dumbbells are predominantly in the plastic crystal phase. Our system¬ 
atic investigation revealed the expected structural changes with larger elongation 
with respect to the hard sphere reference case and very localized orientational cor¬ 
relations, typically just involving next-neighbor correlations. These relatively weak 
correlations are also reflected in only minor changes in the translational and rota¬ 
tional diffusion coefficients for most of the investigated elongations, except for the 
highest elongation in the fluid phase where the rotational diffusion drops by ca. 
30 % when compared to the free rotation. However, the linear response shear vis¬ 
cosity exhibits a dramatic increase even in the plastic crystal phase at high packing 
fractions {4> > 0.5) beyond a critical elongation of about L* = 0.15. This result is 
surprising in view of the relatively weak effects of elongation found before. Here one 
should rationalize that the linear response shear viscosity expresses collective relax¬ 
ation time scales and not (more local) single molecule dynamic properties. Appar¬ 
ently beyond a critical, surprisingly small anisotropy, newly occurring rotational- 
translational couplings must be made responsible for the slow time scales appear¬ 
ing at higher elongations in the crystal US]. Hence, more detailed calculations and 
modeling of these rotational-translational correlations in weakly anisotropic model 
systems shall be interesting for future studies, in particular beyond linear response 
where probably more substantial dynamical effects of increased anisotropy may 
occur p6] . 
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